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ABSTRACT 

The orbital structure of triaxial models with weak central density cusps, 
p oc r~ 7 , 7 < 1, is investigated. The stability of the x- (long-) axis orbit - and 
hence the existence of box orbits - depends sensitively on 7; the range of model 
shapes for which the x-axis orbit is stable becomes progressively smaller as 7 
approaches one. The banana and fish boxlets in the x — z (long axis-short axis) 
plane are stable over a wide range of model parameters. The boxlets in the 
x — y and y — z planes are generally vertically unstable. 



1. Introduction 

The phase space of Stackel potentials is completely occupied by four major families 
of regular orbits: the boxes, and the three families of tubes ( [Kuzmin 1973| ; |de Zeeuwj 



fc hyndcn-Bcll 1985 ). Roughly the same is true in many nonrotating triaxial potentials 



corresponding to mass models with smooth cores: most of the orbits have shapes that can be 
identified with one of the four families of regular orbits in Stackel potentials ( [Schwarzschild 



|1979| ). However the situation is very different in triaxial models with divergent central 
densities, or cusps. While the tube orbits are not strongly affected by cusps, the boxlike 
orbits - orbits with stationary points and filled centers - are often rendered stochastic, 
except when they are associated with a resonant orbit that avoids the center QGerhard fc 



Binney 1985| ; [Lees fc Schwarzschild 199~2| ; |Merritt fc Fridman 1996| ; Papaphilippou fc Laskar 



19971) . 



Kormendy pointed out already in 1985 that many early-type galaxies and bulges have 
central brightness profiles that deviate slightly but systematically from that of an isothermal 
core. The significance of this deviation was not recognized for ten more years, however, due 
to an optical illusion associated with projection onto the plane of the sky. A luminosity 
density profile that varies as r -7 at small radii generates a power-law cusp in projection only 
if 7 > 1. When 7 = 1, the surface brightness exhibits a curving, logarithmically-divergent 
central profile (e.g. Dchnen 1993| , Fig. 1), and for 7 < 1 the central surface brightness is 



finite. The observed brightness profile of a galaxy like M87, which has a cusp with 7 « 0.8 
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( [Lauer et al. 1992]) , differs only subtly from that of a galaxy with an isothermal core 
( [Kormendy 1985| ). This is true in spite of the fact that the cusp in M87 is well resolved 
from the ground. 

Thus even if galaxies were distributed uniformly over 7, their central brightness profiles 
would appear to fall into one of two distinct categories: the "cores" (7 < 1) and the 
"power-laws" (7 > 1). Just such a dichotomy was proposed following the first nuclear 
brightness measurements from HST ( [Ferrarese et al. 1994j ; [Lauer et al. 1995| ). Merritt & 



Fridman (1995) suggested that the centers of all early-type galaxies, including the core 
galaxies, might contain power law cusps in the space density. They used a nonparametric 
deprojection algorithm to confirm this hypothesis for six galaxies observed with HST. 
Gebhardt et al. (1996) subsequently reanalyzed the full sample of 42 early-type galaxies 
observed by Lauer et al. (1995) and verified the power-law nature of the deprojected cusps 
in each case. They also confirmed a tendency, noted earlier by Kormendy (1985), for the 
brightest galaxies to have the shallowest cusps. 

Most of the theoretical work on orbital motion in triaxial potentials has focused on 
mass models with cores, 7 = 0, or with steep cusps, 7 = 2. But galaxies with cores probably 
do not exist; and only the faintest ellipticals, M v > —19, have cusps that are predictably as 
steep as r~ 2 . Here we present the first detailed study of periodic orbits in triaxial models 
with weak cusps, < 7 < 1, characteristic of bright elliptical galaxies. The value 7 = 1 
is a natural one for separating "weak" from "strong" cusps, since a central density that 
increases more rapidly than r -1 implies a divergent central force. For instance, a spherical 
galaxy with Dehnen's (1993) density law 

(3-7)M/rW r\-( 4 -T) 



Ana 3 \aj \ a 

has a gravitational force 



<9<I> GM fr\ 1 -^ ( r V~ 3 

1 + - • (2) 



dr a? \aj \ a 

Dehnen's law, with 7 = 2 flJafi'e 1983Q and 7 = 1 ( [Hernquist 1990|) , has been shown to 
accurately describe the brightness profiles of a number of early-type galaxies. The triaxial 
generalization of Dehnen's law, 

p(m) = ?~ 7 7 ^ ~ 7 (1 + m)- (4 " 7) , (3) 
Anabc 

x 2 iP 1 

m 2 = - + ^ + -, a>b>c>0, (4) 
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was presented by Merritt & Fridman (1996), who derived expressions for the gravitational 
potential, force and force gradients for < 7 < 2. 

Our primary concern here is with the closed "boxlets" , the low-order resonant orbits 
that act as parents of the boxlike orbits in non-integrable triaxial potentials. Most 
important is the x- (long-) axis orbit, which generates the box orbits in Stackel potentials. 
Its instability usually occurs through the bifucation of a 2 : 1 resonant orbit, or "banana." 
(We adopt Miralda-Escude & Schwarzschild's (1989) scheme for labeling the closed boxlets; 
their Figure 4 illustrates the different types.) Higher-order resonances produce the "fish" 
(3 : 2) and "pretzel" (4 : 3) orbits. Boxlets corresponding to each resonance are expected to 
exist in all three of the principal planes of a triaxial model, at least at certain energies and 
for certain values of the model axis ratios. However little is known about the dependence 
of boxlet stability on model parameters. In particular, the vertical stability of the closed 
boxlets has been very little explored. 

In the present paper, we locate and test the stability of the axial orbits, the bananas, 
and the fish as a function of the three parameters that define Dehnen's models: the cusp 
slope 7; the short-to-long axis ratio c/a; and the degree of triaxiality, expressed in terms 
of T = (a 2 — b 2 )/(a 2 — c 2 ). T = corresponds to oblate-axisymmetry, and T = 1 to 
prolate-axisymmetry. For each set of model parameters, we present results as a function of 
energy; thus we explore a four-dimensional parameter space for each of the orbit families. 
A total of approximately 3.5 x 10 4 orbits were investigated. The numerical scheme for 
evaluating the stability of the boxlets is presented in §2. The results for the axial orbits, the 
bananas and the fish are presented in §3, 4 and 5 respectively. Our results are summarized 
in §6. 

There is increasingly strong evidence for dark mass concentrations, possibly 
supermassive black holes, at the centers of many early-type galaxies ( [Kormendy fc| 
Richstone 19951) . We ignore the (possibly substantial) effect that such singularities would 



have on the behavior of boxlike orbits: partly because the universality of central black holes 
has not yet been established; and partly because the behavior of orbits in triaxial galaxies 
with cusps provides an essential first step toward understanding the additional effects of a 
central black hole. 



2. Linear Stability Analysis 



We applied a numerical version of the standard Floquet-Liapunov theory ( [Floquet 1883 
Liapunov 1892| ) to evaluate the stability of the lowest-order periodic orbits in Dehnen's 
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model potential. Our approach was similar to that of earlier studies (e.g. Magncnat 1982 
Pfenniger 1984| ), with one difference. Instead of integrating a set of orbits whose initial 
conditions were offset by small, finite amounts from those of the closed orbit, we integrated 
the linearized equations corresponding to a set of infinitesimal perturbations about the 
closed orbit. Our results are therefore strictly independent of the numerical amplitude of 
the perturbation. 



Let X(£), V(i) be the parametric representation of a closed orbit. The equations of 
motion are 

dX(t) dV(t) 



V(t), 



-V$(X). 



(5) 



dt dt 

A nearby orbit has coordinates X(t) + x(t), V(£) + v(t). For small x and v, the equations 
of motion are 



d{X + x) 
Jt 

d(V + v) 



V 



(6) 



dt 



= -V$(X + x) « -V$(X) + (xV)(V$(X)). 
Subtracting equations (0) from equations (0), we have 

-(xV)(V$(X)), 



"eft 
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The second derivatives of the potential, , are taken along the closed orbit whose period 
is P. Expressions for the $jj are given by Merritt & Fridman (1996). The stability of the 
closed orbit is then determined by the eigenvalues of the monodromy matrix W, where 



x(t + P) = Wx(t) 



(9) 



and x has been redefined to include velocities. We computed the elements of W by 
numerically integrating the six linearized equations of motion (|8|), starting from unit 
perturbations in each of the coordinates. A 7/8th order, variable timestep integrator - the 
routine DOPRI8 of Hairer et al. (1987) - was used. Initial conditions for the non-axial 
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closed orbits were found by an iterative Newton method, and precise orbital periods were 
calculated using the scheme of Henon (1982). 

Two of the eigenvalues of W are guaranteed to be unity; these correspond to 
eigenvectors lying along the closed orbit. The remaining four eigenvalues satisfy Ai = I/A2, 
A 3 = 1/A 4 ; the stability of the closed orbit is determined by the parameters 61 = — (Ai + A 2 ), 
^2 = — (A3 + A4), with > 2 denoting instability flBroucke 1969| ). For orbits that lie 



in a symmetry plane, b\ and 6 2 correspond to perturbations in the orbital plane and 
perpendicular to the orbital plane, or vice versa. For the simplest orbits, the full set of six 
perturbation equations is not required; for instance, the stability of the axial orbits can be 
evaluated using just four independent perturbations. 

Axial orbits presented one complication, however. While the gravitational force — V$ 
is everywhere finite in Dehnen models with 7 < 1, the second derivatives of the potential 
diverge at the origin as <&u cx |xj|~ 7 . To avoid this divergence, analytical expressions were 
derived for the evolution of the perturbations very near the origin, allowing this region to 
be excluded from the numerical integrations. 

The accuracy of the algorithms was checked in a number of ways. The monodromy 
matrix and its eigenvalues were checked to see that they accurately satisfied all of their 
known, exact properties, e.g. detW = 1, Ai = 1/A 2 , etc. Integrations using an algorithm 
based on finite-amplitude perturbations were carried out for some orbits and compared to 
the results of the linearized code. Finally, numerical surfaces of section were constructed for 
a few models and compared to the results from the stability calculations. 

All orbits in a given model were assigned energies from a discrete grid, consisting of 
20 energy values that corresponded to the values of the potential on the x-axis of a set of 
ellipsoidal shells - with the same axis ratios as the density - that divide the model into 21 
sections of equal mass. Thus, shell 1 encloses 1/21 of the total mass, shell 2 encloses 2/21, 
etc. Henceforth units are adopted such that G = M = a = 1. 



3. Axial Orbits 

The long- (x-) axis orbit is stable at all energies in integrable triaxial potentials and 
acts as the parent of the box orbits ( [Kuzmin 1973| ; dc Zeeuw fc Lyndcn-Bcll 1985| ). In more 



general triaxial potentials, like the ones considered here, the x-axis orbit becomes unstable 
at high energies when its frequency of oscillation falls to 1/2 the average oscillation frequency 
of a perturbation in the direction of the short or intermediate axis ( |Mir alda-Escude 75 



Schwarzschild 1989]) . A bifurcation then occurs, with the 2:1a; — z oy x — y banana orbit 
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branching off. A similar bifurcation occurs from the intermediate (y-) axis orbit, producing 
the y — z banana. 

Nearer the center, the y- and z-axis orbits become unstable at the bifurcation points 
of the 1 : 1 loop orbits. In models with a harmonic core, the loop orbits first appear just 
outside of the core flMerritt fc de Zeeuw 198^ ) . The x — y loop is the parent of the short-axis 



tubes, and the y — z loop is the parent of the long-axis tubes. The x — z loop is generally 
unstable and generates no families of regular orbits (|Heiligman fc Schwarzschild 1979] ) . In 



models without cores, the loops and their associated tube orbits may exist all the way into 
the center. 

Figures 1 and 2 display the stability diagrams for the x- and y- axis orbits in Dehnen's 
potential, for three values of T and four values of 7; a total of 4104 orbits are displayed in 
each figure. Four cases are distinguished: stability (•); instability in the direction of the 
longer (+) or shorter (x) of the two remaining axes; and instability in both directions (•). 

The x-axis orbit is stable at low energies in all models that are not too highly elongated. 
As 7 increases above zero, the range of model shapes for which the x-axis orbit is stable to 
the x — z banana bifurcation becomes narrower, eventually including only nearly spherical 
models (Figure 3). For instance, when 7 = 0.9, the x-axis orbit is stable at the half-mass 
radius for c/a > 0.8. Thus, bona-fide box orbits (which require for their existence a stable 
long-axis orbit) will occur at a wide range of energies only in models that are nearly 
spherical, that have very weak cusps, or both. Figure 4 illustrates how the bifurcation of 
the x — z banana orbit takes place at successively lower energies as 7 is increased. 

As c/a is decreased still more, a second 2 : 1 bifurcation from the x-axis orbit sometimes 
occurs in the direction of the intermediate axis, giving rise to the x — y bananas. Following 
this bifurcation, the x-axis orbit is unstable in both directions. As Figure 1 shows, the x — y 
banana bifurcation occurs most readily in nearly prolate models (e.g. T = 0.8) in which the 
y and z axes are nearly equal in length. For nearly oblate models, e.g. T = 0.2, the x and y 
axes are too similar in length for the x — y frequency ratio to ever reach the required value 
of 2 : 1, and the x — y banana does not appear. 

Interestingly, the x-axis orbit can again become stable when c/a is made smaller still 
- typically less than about 0.45. The return to stability coincides with the appearance 
of the x — z "anti-banana" orbit, a 2 : 1 resonant orbit that passes through the center 
( [Kliralda-Escude fc Schwarzschild 1989|) . In nearly oblate models, this return to stability in 



the x — z plane causes the x-axis orbit to become stable in both directions, at least over a 
limited range in energy and in c/a (Figure 1). The axial orbit becomes unstable once again 
at still lower values of c/a, through the appearance of a 3 : 1 resonance in the direction of 
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the z axis (Figure 5). In nearly prolate models, where the x-axis orbit is also unstable to 
perturbations in the x — y plane, the appearance of the x — z anti-banana orbit leaves the 
orbit still unstable to perturbations in the y direction (Figure 1). 

At very low values of c/a, the stability of the x-axis orbit varies in a rapid and complex 
manner with energy and with 7, as resonances of higher and higher order occur. 

The y-axis orbit (Figure 2) is unstable at almost all energies in the full range of models 
investigated here. Instability in the x-direction is expected due to the 1 : 1 resonance that 
generates the x — y loops; these orbits - and their associated family, the short-axis tubes - 
are present almost universally in Dehnen's models. As c/a decreases, the y-axis orbit also 
becomes unstable through the 2 : 1 bifurcation that produces the y — z bananas; the y-axis 
orbit is then unstable in both directions. At still greater model elongations, the y-axis 
orbit returns briefly to stability in the z-direction through the appearance of the y — z 
anti-banana. 

The stability diagram of the z-axis orbit is not given here. This orbit was found to be 
unstable to perturbations in both the x- and y-directions for almost all values of the model 
parameters, due to the 1 : 1 bifurcations that produce the x — z and y — z loops flGoodman 
|fc Schwarzschild 1981|) . Thus, the long-axis tube orbits - which are parented by the y — z 
loop orbits - are almost always present in these models. The only exceptions occur at low 
energies in models with c/a < 0.2 and 7 < 0.5, where the ratio of oscillation frequencies in 
the z- and x-directions is too great for the 1 : 1 resonance to occur. The z-axis orbit is not 
the parent of any banana family. 

In summary: only the x-axis orbit is stable over an appreciable range of model 
parameters. Instability first appears through the x — z banana bifurcation, which occurs 
at progressively larger values of c/a as 7 is increased. The x-axis orbit returns to stability 
for a narrow range of c/a values, typically around 0.2 — 0.4, through the bifurcation of the 
x — z anti-banana. 

4. Banana Orbits 

Banana orbits first appear as 2 : 1 bifurcations from the x-axis orbit [x — z and x — y 
bananas) or the y-axis orbit (y — z banana). Stability diagrams for the x — z and x — y 
bananas are shown in Figures 6 and 7. Once these orbits bifurcate from the axial orbits, 
they continue to exist (though not always stably) at all higher energies and all lower values 
of c/a, for any given value of 7 and T. 
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The x — z banana exists and is stable over the widest range of model parameters (Figure 
6). It becomes unstable only in nearly-prolate models, through a vertical 2 : 1 bifurcation; 
for example, when T = 0.8, instability first occurs when c/a drops below ~ 0.35. Thus, in 
virtually all models with axis ratios comparable to those of real elliptical galaxies, either 
the x-axis orbit or the x — z banana orbit is stable and can act as the parent of a family of 
regular boxlets. 

The x — y banana exists only in relatively prolate and elongated models, in which the 
x — y frequency ratio can exceed 2 : 1 (Figure 7). For moderate values of c/a, however, the 
x — y banana is generally vertically unstable; its instability follows closely that of the x— axis 
orbit in the z-direction (Figure 1). For smaller values of c/a, e.g. c/a < 0.5 for T = 0.8, the 
x — y banana returns to stability via a 2 : 1 bifurcation at approximately the same model 
parameters where the x— axis orbit also becomes stable in the z— direction (Figure 8). 

The y — z banana was found to be almost universally unstable in the vertical (x) 
direction, mimicking the x- instability of the y-axis orbit. 

The regular orbits associated with the x — z bananas are known to be heavily populated 
in self-consistent triaxial models with cusps (e.g. [Bchwarzschild 1993 ; |Merritt fc Fridm"an| 



1996|) . However their usefulness is limited by the fact that the banana orbits become 
more strongly curved as c/a decreases, so that in sufficiently elongated models, the boxlets 
associated with the x — z bananas are all rounder than the model figure. Figure 9 shows 
the value of c/a at which the x — z banana has an elongation - defined as \z max /x max \, the 
ratio of its extent in the z and x directions - that just matches that of the model. Figure 9 
suggests that the boxlets associated with the x — z banana will be useful in reconstructing 
triaxial models for c/a > 0.4. 

In summary: the x — z banana exists and is stable for the widest range of model 
parameters. It is slender enough to match the shape of the model when c/a > 0.4. The 
x — y and y — z bananas either do not exist, or are vertically unstable, in all models with 
moderate values of c/a; the x — y banana returns to stability only in highly elongated, 
nearly prolate models. 



5. Fish Orbits 

Fish orbits first appear as 3 : 2 bifurcations from the x-axis orbit [x — z and x — y fish) 
or the y-axis orbit (y — z fish). Stability diagrams for the x — z and x — y fish are shown in 
Figures 10 and 11. Once these orbits bifurcate from the axial orbits, they continue to exist 
at all higher energies and all lower values of c/a, for any given value of 7 and T. 



- 9- 



The x — z fish exists and is stable in models with moderate values of c/a (Figure 10). 
In nearly oblate models, the x — z fish first becomes unstable to perturbations in the orbital 
plane, while for strongly triaxial and prolate models, instability first appears in the vertical 
(y) direction (Figure 12). However this orbit exists stably over a substantially narrower 
range of model parameters than the x — z banana. 

The x — y fish is only important in strongly prolate models (Figure 11). In strongly 
triaxial and oblate models, it either does not exist, or is generally unstable to vertical (z) 
perturbations. 

The y — z fish was found to be unstable to vertical perturbations for almost all values 
of the model parameters. 

In summary: the x — z fish is important in moderately flattened, oblate/triaxial models, 
while the the x — y fish is important in highly prolate models. In strongly triaxial models, 
instability in both families first occurs in the vertical direction. Neither orbit exists stably 
for as wide a range of model parameters as the x — z banana. 



6. Discussion 

Our results highlight the sensitivity of boxlike orbits to small changes in the central 
density structure of a triaxial galaxy. The range of model shapes for which the x-axis orbit 
- the progenitor of the box orbits - is stable becomes rapidly narrower as the cusp slope 7 
approaches one. For 7 > 0.8, only triaxial models with moderate axis ratios, c/a > 0.75, 
contain stable long-axis orbits at most energies. Regular box orbits may nevertheless be 
present in the majority of bright elliptical galaxies, since the distribution of intrinsic shapes 



of these galaxies is strongly peaked at c/a « 0.8 (|Tremblay fc Merritt 1996|) , and many 



ellipticals have cusp slopes 7 < 0.8 ([Merritt fc Fridman 1995| ; |Gebhardt et al. 1996|) . 



Of the remaining families of boxlets investigated here, those lying in the x — z plane 
were found to be most important, in the sense that they exist and are stable in models 
with the widest variety of cusp slopes and axis ratios. In particular, the x — z bananas are 
unstable only in highly flattened, nearly prolate models, and they are often more elongated 
than the model figure, making them useful building-blocks for a self-consistent galaxy. The 
boxlets in the x — y and y — z planes are almost always vertically unstable in strongly 
triaxial models (when they exist), and it seems unlikely that they could play an important 
role in maintaining triaxiality in real galaxies. 



We have investigated only motion in the principal planes. Other resonant orbits - 
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although generally of a higher order - exist outside of the principal planes. Merritt & 
Fridman (1996) cataloged the orbits in a single triaxial model with Dehnen's density law, 
7 = 1, c/a = 0.5 and T = 0.5. Their Figure 4 confirms that only the x — z banana and 
the x — z fish, among the closed boxlets investigated here, are generators of regular orbit 
families. The x — y pretzel was also found to be important at low energies, as were a number 
of higher-order resonances outside of the principal planes, including the 4 : 5 : 7, 5 : 6 : 8, 
and 6:7:9 closed boxlets. 

Although instability does not imply stochasticity, unstable resonant orbits are often 
associated with chaos. We expect that the phase space surrounding many of the unstable 
boxlets is stochastic. However, the timescale over which the stochasticity manifests itself in 
the orbital motion is likely to be a strong function of the model parameters, particularly 
the cusp slope 7 and the energy. We will return to this question in a future study. 

This work was supported by NSF grant AST 90-16515 and by NASA grant NAG 
5-2803. 
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Fig. I.- 

Stability diagram for the x-axis orbit. • stable; + unstable in the direction of the y 
axis; x unstable in the direction of the z axis; • unstable in both directions. 

Fig. 2.- 

Stability diagram for the y-axis orbit. + unstable in the direction of the x axis; x 
unstable in the direction of the z axis; • unstable in both directions. 

Fig. 3.- 

Stability region for the x-axis orbit. Solid lines correspond to models in which the 
x-axis orbit is just unstable, at the half-mass radius (shell 10), to bifurcation of the x — z 
banana orbit. Box orbits do not exist outside the half-mass radius in models below these 
lines. 

Fig. 4.— 

Stability characteristics of the x-axis orbit when perturbed in the z-direction, for 
c/a = 0.6 and T = 0.5. Frequencies and e-folding rates are expressed in units of the orbital 
frequency. Bifurcation of the x — z banana orbit occurs at successively lower energies as the 
cusp slope 7 is increased. 

Fig. 5.- 

Stability characteristics of the x-axis orbit when perturbed in the ^-direction, for 
7 = 0.4 and T = 0.5. Frequencies and e-folding rates are expressed in units of the orbital 
frequency. The x-axis orbit returns to stability through bifurcation of the x — z anti-banana 
orbit. 

Fig. 6.— 

Stability diagram for the x — z banana orbit. • stable; x vertically unstable. No 
symbol indicates that the x — z banana does not exist. 

Fig. 7.- 

Stability diagram for the x — y banana orbit. • stable; x vertically unstable. No 
symbol indicates that the x — y banana does not exist. 

Fig. 8.- 

Stability characteristics of the x — y banana orbit when perturbed in the z-direction, 
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for 7 = 0.6 and T = 0.8. Frequencies and e-folding rates are expressed in units of the orbital 
frequency. 

Fig. 9- 

Value of c/a at which the elongation of the x — z banana matches that of the model, 
for T = 0.5, and for four values of the cusp slope: 7 = 0.2, 0.4, 0.6, 0.8. 

Fig. In- 
stability diagram for the x — z fish orbit. • stable; + unstable in the plane; x vertically 
unstable. No symbol indicates that the x — z fish does not exist. 

Fig. 11.— 

Stability diagram for the x — y fish orbit. • stable; + unstable in the plane; x vertically 
unstable. No symbol indicates that the x — y fish does not exist. 

Fig. 12.- 

Stability characteristics of the x — z fish orbit when perturbed in the ^/-direction, for 
7 = 0.8 and T = 0.8. Frequencies and e-folding rates are expressed in units of the orbital 
frequency. 
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